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Abstract 
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Flexible boundary condition methods couple an isolated defect to a harmonically responding 
medium through the bulk lattice Green's function; in the case of an interface, interfacial lattice 
c/3 \ Green's functions. We present a method to compute the lattice Green's function for a planar in- 

terface with arbitrary atomic interactions suited for the study of line defect/interface interactions. 
The interface is coupled to two different semi- infinite bulk regions, and the Green's function for 
interface-interface, bulk-interface and bulk-bulk interactions are computed individually. The elas- 
tic bicrystal Green's function and the bulk lattice Green's function give the interaction between 
bulk regions. We make use of partial Fourier transforms to treat in-plane periodicity. Direct in- 
version of the force constant matrix in the partial Fourier space provides the interface terms. The 
general method makes no assumptions about the atomic interactions or crystal orientations. We 
simulate a screw dislocation interacting with a (1012) twin boundary in Ti using flexible bound- 
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conditions give the correct core structure with significantly less atoms required to relax by en- 
ergy minimization. This highlights the applicability of flexible boundary conditions methods to 
modeling defect /interface interactions by ah initio methods. 
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I. INTRODUCTION 

Accurate atomic scale studies of lattice defect geometry is the key to any modeling of 
their effects on material properties. However, the long-range (elastic) displacement field of 
isolated defects, e.g., dislocations, is incompatible with periodic boundary conditions typi- 
cally used in computer atomistic simulations. Fixed boundary conditions require simulation 
sizes large enough for the elastic solution to be accurate — a size typically beyond even mod- 
ern density-functional theory methods. Flexible boundary condition methods avoid these 
issues by relaxing the atoms away from the defect core through lattice Green's function 
(LGF) as if they are embedded in an infinite harmonic medium. Hence, the atomic scale 
geometry of the defect core is coupled to the long-range strain field in the surrounding 
medium. Sinclair et al. introduced flexible boundary conditions for studying defects in bulk 
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assical potentials and 



materials[l|] such as cracks|2|, |3|, dislocations [#i7|, vacancies with c 

isolated screw or edge dislocations with density-functional theory |8l-lll|. Flexible bound- 
ary conditions use the LGF corresponding to the specific geometry of the problem. For 
instance, line defects in the presence of interfaces require the interfacial lattice Green's 
function (ILGF). Line defects in interfaces affect the mechanical properties of composites. 



two-phase or polycrystalline materia 



s where heterophase or homophase interfaces interact 



with defects. Tewary and Thomson 12| proposed a Dyson-equation calculation of the in- 
terfacial lattice Green's function suitable for materials with short-range atomic interactions 
and simple crystal structures. We present a general — for all types of interactions and in- 
terface orientations — accurate method to compute the interfacial lattice Green's function, 
suited to use in density functional theory. Specifically, this method is applicable to studies 
of line defects interactions with planar interfaces such as disconnections in interfaces and 
dislocation or crack tips interacting with grain boundaries and two-phase interfaces. We 
compute the Green's function for a (1012) twin boundary in Ti to simulate a screw dislo- 
cation interacting with the twin boundary using flexible boundary conditions. Section HTl 
reviews the harmonic response functions: the force constant matrix and the lattice Green's 
function. Section IIIII explains the general procedure for evaluation of the interfacial lattice 
Green's function and section IIVI applies the method to modeling the interaction of a screw 
dislocation with Ti (10l2) twin boundary. The end result is a computationally tractable, 
general approach usable for studies of defects in interfaces. 



II. HARMONIC RESPONSE 



Harmonic response is characterized by a linear relationship between forces and 



displacements 13|] . Lattice Green's function G{R, R') relates the displacement u{R) of atom 



R to the internal forces f{R') on another atom R' of the crystal through 

t7(^) = -X^G(^,^')/(^')- (1) 

R' 

Conversely, the forces on an atom can be expressed in terms of displacements through the 
force constant matrix D.{R, R') by 

/(^) = -X^D(i?,i?')^(^')- (2) 

R' 

Translational invariance of an infinite crystal makes G and D. functions of the relative posi- 
tions of the atoms. Substituting Eqn. (ED into Eqn. (^ gives ^G{R - R')D{R') = 16 (R), 

R' 

where S{R) is the Kronecker delta function. A constant shift in atom positions does not 

produce internal forces; hence, J2rI2.{R) = 0, and so G{R) is the pseudo inverse of D{R) 

in the subspace without uniform displacements or forces. In a bulk geometry, the Fourier 

transform of the lattice functions are defined as 

,,__ ^ ^ ^ ^ r (j'iu 

G{k) = Y,e^'-^G{R), G{R)= 7^^e-^^-«G(A;) 

^ Jbz {^T^y 

R 

where the summation is over lattice points. In reciprocal space, the matrix inverse rela- 
tion G{k)D{k) = 1 and the sum rule D{0) = require that G{k) has a pole at the F-point. 
While computation of the force constant matrix D(i?) — and subsequently D{k) — is straight- 
forward, G{R) can not be computed directly due to its long range behavior. Instead, we 
invert D(k) to get G(k) and then perform an inverse Fourier transform. Convergence of the 
inverse Fourier transform requires an analytical treatment of the pole at the F-point |14J.|15|. 
In an interface geometry, translational invariance is broken in the direction perpendicular 
to the interface; we use Fourier transforms in the interface plane only. This produces an 
infinite dimensional dynamical matrix that can not be simply inverted, but requires a more 
complex computational approach. 



III. COMPUTATION OF LATTICE GREEN'S FUNCTION FOR A PLANAR IN- 
TERFACE 

Figure [T^ shows two lattices, A and /i joined at a planar interface. Each set of vectors 
d[ '^, d^ '^ and 03 '^ give the periodic directions in their corresponding lattice. We introduce 
integer matrices M and M ^ and deformation operators F_ and F_'^ so that 

F_'^[ai'^,a2'^,a2,'^\K = Ti ,T2 ,T3 (3) 

^A -» fi -t -»A -* li -t 

to define the supercell. We use Ti = Ti = ti and T2 = T2 = ^2 as nonparallel vectors to 
define the interface plane where ^2 will be the periodic threading vector for a line defect in 
the interface. The combined lattice has translational invariance in ti and ^2 directions in the 
interface plane while the periodicity is broken in directions outside the plane. Introducing 
a threading direction reduces the problem to 2D (i.e plain strain or anti-plane strain condi- 
tions). We confine our calculations to the plane orthogonal to ^2 and define the Cartesian 
coordinate x, y, z so that ti-x = Oq, t2 = |t2| y and z = x xy. Note that in general Qq ¥" |^i| 
because ti and ^2 can be nonorthogonal. Specifically, the lattice positions, R = xx + zz and 
the Fourier vectors, k = kxX + kzZ, will be 2D vectors through out this paper and 

D^^,{R,R') = D^A^-x'-z,z') 

with a and a' identifying the xyz components of the second rank tensor D in Cartesian 
coordinates. We index atoms in our computational cell with integer I at position {xi,zi)] 
due to periodicity in the x direction, each atom also occurs at Xi + naox for integer values 
of n. The partial Fourier transform is 



n=— 00 

D^A^i - ^v + na,- zu z,) = ^ e-^'=^(^'-^"+"'^°)D,;,„,z,(A;,)cifc, (4) 

for all pairs /, /'. Note that "/" indexes layers of atoms with particular z values. There may 
be two different layers that have equal z: zi = zi' while / 7^ /'. D{kx) is infinite dimensional 
due to infinite values of /. 

To avoid the inversion of infinite dimensional D{kx), the geometry is divided into two 
semi-infinite bulk regions coupled with an interface region. Figure [T]d shows the schematic 
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FIG. 1: (a) Bicrystal /i and A, (b) separation into bulk and interface regions and (c) the Ti (1012) 
twin boundary . Two different lattices, A and /x are connected through a planar interface. The unit 
cells of A and ^ are given by di '^,02 '^ and 03 '^ — all of which must be lattice vectors in A and 
/i. The combined lattice has the periodicity of the interface in ti and ^2 directions. Introducing 
a line defect threading direction t2 reduces the problem to 2D in the plane normal to t2- In (b), 
the crystal is divided into two semi-infinite bulk regions, bulk A and bulk /x symbolized by (+) and 
(— ) respectively, coupled with an interface region (I). The bulk regions are far from and affected 
only through an elastic effect by the interface. The force constant matrix between atom pairs in 
the bulk is not affected by the interface. The remaining layers are included in (I), (c) shows the 



periodicity vectors for the Ti (1012) twin boundary. The interface is defined by ti = VSa^ + c^x 
and t2 = CLy where a and c are the hep unit cell parameters in Ti for both A and fi. /i is the 
reflection of A about the interface plane. 



divisions of the regions in an interface geometry consisting of lattices A and /i. The "bulk" 
regions represent layers of atoms that are far from and affected only through an elastic field 
by the interface. The atomic scale interaction between atom pairs are as if they were in 
their corresponding bulk geometry. Bulk A and bulk /i are symbolized by (+) and (— ) in 
our notation. The remaining layers, affected by the reconstructions near the interface, are 



included in the "interface" region (I). We define the interface region as atoms where the force 
constant matrix differ from those in the bulk lattice. For specific geometries, additional bulk 
layers may be included in the interface to insure a smooth transition between the regions. 
We block partition the infinite dimensional Dai^a'i'{kx) and Gai,a'i'{kx) based on the atom 
region (+, — , or I) of indices / as 



Dih) 



( D'\kx) 


D'-{K) U'+{kx)\ 







(5) 



where / > /+ belong to (+) region, / < /_ belong to (— ) region and the finite-dimensional 
region is (I). D{kx) and G{kx) are Hermitian and satisfy 

2_^ 6al,a"l"{kx)Ga"l",a'l'{kx) =Saa'Sll'- (6) 



a"l" 



We construct D{kx) by direct calculation of D^^iixi — x'i+ nao; zi, z[) followed by a partial 
Fourier transform according to Eqn. (|1]) and block partitioning as in Eqn. ([5]). Note that due 
to the finite number of interface layers and decay of the force constant matrix, the infinite 

dimensional non-zero sections of D{kx) consists of , — h and ++ interactions (bulk-like 

regions with themselves) which we explicitly avoid in our approach. 

The infinite dimensional blocks of G{kx) are known from bicrystal elastic and bulk lattice 
calculations. The distance between + and — is large enough for the elastic Green's function 
to be applicable; the real space solution of G^^ is calculated from the bicrystal elastic Green's 
function in both plane strain and anti-plain conditions proposed by Tewary et. al 16|. We 
partially Fourier transform the real space solution by a continuum version of Eqn. (jl]). 



GaLa'l'\kx) 



G-Z{x-Zuzv)e"'^''dx. 



-qq' V"^' ^'' 



(7) 



G^^ is the conjugate transpose of G^^ due to G{kx) being Hermitian. The functional form 
oi G^^{x; zi,zi') consists oi real parts oi\n{x+PgZi+p'^,Zi') where Pg and p^, are the complex 
roots of the sextic equation of anisotropic elasticity for bicrystal A/i and g, g' = 1, 2 in plain 
strain and 1 in anti-plane conditions IGj. We rewrite ln(x + 'jfl, + iPf^,) with 



The Green's function in real space is the real part of the complex logarithm with the form 



G;+(a:; z^ zy) = J^A In {x + il^f + (A^,') 



+ &11' arctan 



pqq 
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,^ + rCi' 



where a^'^, and 6^^, are real valued coefficients of the term qq' . Eqn. ([8]) is obtained by 



rewriting Eqn. (60) in 16|. The partial Fourier transform is 



G-+ (k ) - ^ V 2a'"'' e"l''"''4-'^"''^ + tb'"'' -^e'l'''''^-'^"''^ (9) 

I ^ I I J. I qqi I I 

with a first order pole at k^ = 0. The 1/ |t2| prefactor is required for the elastic and lattice 
Green's functions to have consistent units of (lengths/energy). We separate the pole from 
the remainder of the Green's function 

Q- + 

Gata'Af''^) = iff +G:!^'Ak.). (10) 

The pole with a constant coefficient G , = — -rt T^o „/ o-ll/ will be treated analytically while 
the nonsingular remainder G ,,,(^x); will be treated numerically. 

The G (k^) and G^^^k^) blocks in Eqn. (EJ are obtained from the bulk lattice Green's 
function of A and /i lattices plus an elastic term due to the presence of the interface. The 
full Fourier transform of the bulk LGF G'^'^{k) is the inverse of the bulk dynamical matrix 
from Section [TTl The partial inverse Fourier transform gives the Green's function in terms 
of k^ and atom indices 



G7i.a'v{K) = ^ r ^ GZ>{k)e-''^^^^-'"^dk, (11) 



for k = {kxX, kzZ) in the Brillouin zone (BZ), Abz the area of the BZ and ki{kx) and kf{kx) 
showing the initial and final values of kz at each k^- G'Z,{k) has a second order pole at 
k = \/k'^ + k1 = which is responsible for the logarithmic long range behavior of LGF in 
real space. The LGF in reciprocal space is 






GZ'ik) = '^^Uk) + GZik) 

where G is the k direction-dependent elastic Green's function and fc{k) is a cutoff function 
that vanishes smoothly at the edges of the BZ. In general anisotropic cases, G (k) is 

represented by a Fourier series expansion as G ,{k) = >^ G \ e*"'^'= where 0^ is the 

n=0 

angle of k relative to an arbitrary in-plane direction and the truncation A^max is sufficiently 
large [l4|. The integrand in Eqn. flTTl) is not singular for fc^ 7^ however the k'^ pole in G°'°'{k) 



results in a pole of order \kx\ in Q^cda'vQ^x)- To treat the small k^ behavior analytically, we 
integrate Eqn. flTTj) as four terms 

fkf(kx) rkf{kx) 

Ikiikx) Jki{kx) 



/ GZ'{k)e-'^^^''-'''^dk, = / GZ'{k){e-'''^^"-'''^ - l)dk, (12) 

J kiikrr) Jki{kx) 

+ / ' GZ.,{k)-^^^,f.{k.,k,)dk, 

Jkdkrr.) '^X ' '^Z 



+ / T^fT^ifcik.,k,)-l)dk, 

Jkdkx) ki + ki 

Jk,(kx) ki + ki 



where G ' is the n = coefficient in Fourier expansion of G ,{k). The first three terms 
in Eqn. ( TT2l) are evaluated numerically while the last integral is 

r^ ^^ G^o,! ,, _ '^Qa,a> ^--,0 / arctan(fc/(/c^)//c^) - arctan(fc^(/c^)/A;^) _ tt \ 

/ , , /p2 I /p2'^'^^ ~ l/r I ~aa' I k \k \ I 

Jkiika:) '^x^'^z \'^x\ \ r^x \'^x\J 

(13) 
where i^"? is the pole and the remaining terms are added to the numerically evaluated 
part. We add an elastic correction term to QZ'ikx), due to the interface obtained from Eqn. 
(59) in [l6|. Combining Eqn. (TT3|) . Eqn. ( !T2|l . and Eqn. (TTOj) produces 



..era 



G 

G7Ukx) = =rff^ + G:,^,{k.)- (14) 



^x\ 



Eqn. (J5]) has unknown blocks G^^{kx), G^^{kx). Direct substitution of the block partitions 
gives 



G'^ik,) = -{D'\kx))-' J2 D'"\k.x)G"'"{kx) (15) 

G'\k,) = {D'\k,))-'+ Y.^G'\kx))-'D''^{kx)G''''\kx)D'^'\kx){D'\kx)r\ (16) 

u'a=± 

Note that by choosing the appropriate set of independent equations we manage to avoid the 
calculation of the infinite dimensional D"^" (kx)- The finite range of D^'^ikx) means that only 
a finite subset of atoms in each semi-infinite ± region are considered for G"" (kx). To treat 
the poles in G^^{kx) and G^'^{kx) analytically, we use a kx expansion oi D{kx) = D + D{kx) 



derived from Eqn. (jlj) where D{kx) = D^k^ + 0{kl) . Therefore, for small kx 

{Dik.yr' = [5+5(m]"' (17) 

= jy^ - k^D'^D^D'^ + 0{kl). 

Using the small k^ expansions for the bulk Green's functions with Eqn. f lTSjl and Eqn. f lTB]) 
gives 

G'''{kx) = -^g'" + G'\kx) and G'\kx) = -^g'' + G'\kx) (18) 

where 

G =-{D ) 2IQ G andG = {D )~^ + 2^(5 ^Q G D {D Y^ 

are the constant coefficients of the pole and G {kx) and G (kx) include the remaining 
nonsingular terms. G (kx) and G {kx) have a cusp approaching kx = Q and the value at 
fca; = is 



g'\q) = -ib'')-' Y, ti'' G'''\^) (19) 

(t'=± 

G'\<d) = {d")-'+ J2 (5'V'5"^G'^'^'(0)5^''(5'V' (20) 

a,cr'=± 

where G (0) is calculated in Appendix |Al To ensure a smooth transition between interface 
and bulk regions, we compare the pole terms and the cusps for atom indices at the boundary 
between the regions (i.e /+ and /_). Labeling /o- as cr = ± or (I) does not change the material 
response. Specifically we should have 

g'] „ = G''\ , g'] „ (0) = G7 „ (0) (21) 

Eqn. f l?I]) determines the finite size effect in the interface. Note that once the bulk force con- 
stant matrix is known, identifying atoms in the interface region does not require additional 
computation effort. 

Evaluating the Green's function in real space between to atoms (x/ -|-nao, zi) and [xv, zy) 
requires a partial inverse Fourier transform over Eqn. ( !T8|) . 

G'Sa'/'(^x)e-''=^("'-""+""")rfA;, (22) 

"fvmax 



and 



G^a'i^i -^i' + nao;zi,zi, 



GZc.'i'{k.y''^^^'"''^''"°^dk^ 



(23) 



The G term in Eqn. ( ITSll is treated analytically via 

r"00 -| 



% 



e ''"'"''dkr 



-2 In \x\ . 



Therefore 



e '"'''' dk^ = -2 In |a;| + 2Ci{k^i,^x). 



(24) 



Note that lim —2 In |a;| + 2Ci = 27 + 21n(A;jnax) where 7 ~ 0.577215 is the Euler constant. 
The partial inverse Fourier transform for G terms are evaluated numerically over a discrete 
kx mesh of size N^ 



G[l,i^i-xi'+nao;zi,z,) = j^J2G'j^^^,^Xk.)e ^'^^^^ "" 



+nao) 



and 



G'^Jxi - XV + nao; z^ zv) = jr-H Gai,a'i'(^-)^ 



-ikx(xi—Xii+nao) 



(25) 



(26) 



Table [T] summarizes the method. 



IV. APPLICATION: LATTICE GREEN'S FUNCTION FOR Ti (1012) TWIN 
BOUNDARY 



We use the method to compute the ILGF for a Ti lattice containing (1012) twin boundary. 
The geometry of this boundary is shown in Figure [It. The F_ ''^ and M '^ matrices are 



F^'/^ = I, M^ 



2 
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1 





1 





1 
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2 
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1 





-1 





1 



The twin boundary is defined by ti = v3a^ + c?a; and ^2 = o,y where a and c are the hep 
lattice constants in Ti. Lattice fi is the refiection of A about the twin boundary plane. The 
force- const ant matrices D{R) are computed using LAMMPS package jl7| with a Ti MEAM 
potential with the maximum cut off distance of 5.5A 18|. The partial FT in Eqn. (jl]) is done 
by a uniform discrete mesh of 40 k^ points over (— 7r/ao, vr/ao) where Oq is the periodicity of 
the geometry in x direction and equal to |ti| in this case. The same kx values must be used 



10 



TABLE I: Summary of the procedure for ILGF computation. Regions (+, — , and I) are defined 
in Figure [Hj. G ^{x;z,z') is the elastic Green's function for a bicrystal computed by Tewary 
et al.ll(y\. G'^'^ik) is the LGF in bulk a = it. FT prefactors required to maintain the consistency 
between elastic bicrystal GF and bulk LGF solutions are also listed. 

1. Compute D_aa'i^i ~ ^i' + ^^o; ^u ^u) directly. Divide the geometry into — , I, + regions. 

oo 

2- Bai,c.'v{kx) = Y. e'''^^'"~''''^''''°'^Dil,{xi-xi>+nao;zi,zi,), ao = periodicity in x direction 
and cr = ±,L Eqn. (gD 

/oo 
G~+{x;zi,z,)e''^-dx, h = ^, 61&2 = Abz- Eqn. © 
-00 

Jki(kx) 

5- GS5;1,,(A:.) = :=gf + GZ^^,,{k.). Eqn. d 



fi'^ , ^l<. 



(HD 



'^- Gi'ia'l'i^ = ^l-^l' +'^O'0;Zi,Zi') 



-'^(-21n|x|+2Ci(^x))+^ ^ Gll^,,i^)e-'^\ Eqn. m-m 



m=l 



in (+),(—) and (I) regions. Limits of k^ in Eqn. ( ITT]) are then chosen so that the equivalent 
of Abz is covered in both (+) and (— ). The first three integrals in Eqn. flT^ are evaluated 
numerically over a uniform kz mesh of 160 points at each k^. For l/c^l < 0.l7r/ao, the density 
of kz mesh is doubled to insure the convergence around the discontinuity at F-point [M, [l5| . 
Figure [2] shows the supercell with bulk (+/— ) and interface (I) divisions and the paths 
along which LGF is evaluated for testing purposes. G^xi^ — xi']Zi,zi') is plotted along a 
vertical and six horizontal paths in the supercell where the reference atom /' is the first atom 
{xii = 0) in the horizontal paths and the atom right below the interface in the vertical path. 
Bulk response along zi = z^ paths is gradually recovered as the paths get farther from the 
interface and closer to the (— ) region. In addition, it is worth noting that paths 1 and 2 
are located in bulk and interface regions respectively. Therefore, the LGF is obtained from 

11 



' along horizental paths 




FIG. 2: [1210] projection of the Ti supercell containing a (1012) twin boundary. The supercell is 
divided into bulk (+/— ) and interface (I) regions, y axis is pointing into the plane. Variation of 
the Gxx component of the lattice Green's function is plotted along six horizontal and one vertical 
paths. The reference atom {x',z') is the first atom in horizontal paths and the atom right below 
the interface in the vertical path. Bulk behavior along the z = z' paths is recovered away from the 
interface. The long range behavior of the LGF matches the EGF along the vertical path, while 
deviating for small z — z' . 



the bulk lattice Green's function along path 1 and from the ILGF method along path 2. 
The good agreement between the response of these two paths verifies the smooth transition 
between the bulk- interface divisions. G^xi^ — x'; z, z') as a function of z is also plotted for 
atoms along the vertical line shown on the supercell in Figure El The reference atom is 
located on the vertical line at xu = x', zu = z' = — 1.413A which is right below the interface. 
The long range behavior of the ILGF matches the EGF. 

We apply the computed ILGF to simulate the interaction of a [1210] screw dislocation 
with the Ti(1012) twin boundary by flexible boundary conditions [l 



nj with a Ti MEAM 



12 



potential [18i]. Periodic boundary conditions are applied along the dislocation line. Flexible 
boundary conditions relax atoms surrounding the dislocation core region with the lattice 
Green's function as if they are embedded in an infinite medium. Conjugate- gradient method 
relaxes the atoms around the dislocation core (region 1). This process generates forces 
on atoms of the intermediate region (region 2). ILGF relaxes the forces on region 2 and 
updates the positions of the outermost atoms (region 3), originally obtained from the elastic 
displacement field of the screw dislocation. To verify the results, we also modeled the same 
dislocation/interface geometry with fixed boundary conditions using supercell radii of 12- 
50b; b is the magnitude of the Burgers vector equal to 1^21 • Outer layers of atoms in a region 
of width 3b are frozen to elastic displacement field of the screw dislocation and the inner 
atoms are relaxed through the conjugate-gradient method using Ti MEAM. Large supercells 
are required to minimize the effect of free surfaces created by the fixed boundaries. 



Figure [3] shows the differential displacement maps[19|] of the screw dislocation core struc- 
ture in the Ti (1012) twin boundary obtained by fixed and flexible boundary. Fixed boundary 
conditions result in a finite size effect that is removed with flexible boundary conditions, or 
with significantly larger calculations. For supercell radii i? < 17b (i?=17b corresponds to 
1312 atoms relaxed), the dislocation center is trapped in the interface while for R between 
18 and 50b — corresponding to 1474 and 11364 atoms respectively — the dislocation center 
moves out of the interface towards the bottom lattice. This is possible due to the broken 
mirror symmetry at the twin boundary for this MEAM potential. The flexible boundary 
conditions supercell has R=12h and 652 ((1):73, (2):219, (3): 360) atoms. The core structure 
from flexible boundary conditions is in good agreement with large fixed boundary conditions 
results- hence the correct structure can be obtained using flexible boundary conditions with 
significantly less atoms than with fixed boundary conditions. 

V. CONCLUSIONS 

We developed an automated computational approach to calculate the lattice Green's 
function of crystals containing planar interfaces for arbitrary force constants and interface 
orientations. This method is more general than the previous Dyson-equation approaches in 
the sense that it can consider long range atomic interactions and reconstructions near the 
interface. We computed the ILGF for a Ti (1012) twin boundary with a Ti MEAM potential 
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Fixed BC Fixed BC Flexible BC 

(R=I2-17b) (R=I8-50b) ((l):5b, (2):4b, (3):3b) 

04^ o 435 o 43; 

' 3,% • iii.% ' 3,% • ' 2,% • as ' 2,% • ' 2,% • as 3,% • 

5/.* - ^'-* 5v% 4,%* 4-* ^ 5,% 5^^» . ^3-* 5,% 3^^* 4,% ^ • 4,% s;,* . ^3-* 5,% 4^^* 4,% ^ • 4,% 

6>°4.% lA%6>°5,% 3-%\% 03,% 6>°3.% ^^ " •* 5> ° S,% ^-^ \% ,03.% 7> ° 4,% ^z* ^^ 5> ° 5,% 2.%\% 03^^ 

, 9(% • lis , 6.% • as , 5,% • , Sffe • LZS , 7ffe • 7S , 4,% • , 81% • 12S , 7ft • IS , ' 5,% 1 

2.% ' 4ft 2.%ll,% H* ■ 10% 1% 5.% 2.% • 2.% ' 5,% 2.% s/h* 2.% „^ ^^» 4.% 2.% » 2.% ' 5> 2,%l(^% R* 2.% g^. ^^» 4.% 2.% 

o 'o*'o o '00 O '00 

3,% 2.% o 3^^ 2^„i^% 6_^ o ^^ , 3.% 3,% 3.% o 3_^ 2^„i^% 6_^ o ^^ , 3,% 

^ -T- • "- - • ■ • 1|* Z37% « 9(% , 2.% » , • • 1}* Z3B%-» 9fi , a.% • 

'%17% 1-* ^ 6^ \% 2.% \% 5^ 4.% UjJftl^T^ 3,% I 6> \^ 2,% \% 5^ 5.% 1 l^.'^T^^ 3,% ' 6> 'g,^ 2.% 

V / O ^ \ O Q \_ O 

2.* 3% o ^„ as, o .„ 4.% H* 1 ijk<» ,„ 5Si 3.ft a!k o , 




4iS 4% ■ \ '-^ 1% ^T ^'° 3,% 5,ii 2,% 3.° V liJft^3Si <>/« ^ 3,1; 2.% 'v* 4.% *'° \ 1>« 2.% ^/^ 4.-« 2.% 

• 7ft . l£a- • 5ft . lift • 8ft . • 9ft ,3AV» 11% . 4.% • 6ft . • 9ft .3iliJ» ia% _ 2.% • 6ft 
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o^4So o^i%o o 5fto 



FIG. 3: Differential displacement maps of a screw dislocation core in Ti (1012) twin boundary 
computed by fixed and flexible boundary conditions. Fixed boundary conditions cause a supercell 
size effect which is evident from different core structures for radius R smaller or larger than 17b 
(1312 atoms relaxed). Flexible boundary conditions give the same core structure as the large fixed 
boundary conditions supercell with significantly less atoms required to relax by energy minimization 
(i.e 73 atoms in region (1) and 652 atoms total). 

and studied the screw dislocation/twin boundary interaction using flexible boundary con- 
ditions. Our results show that the ILGF flexible boundary conditions method predicts the 
correct dislocation core structure. Moreover, the energy minimization stage of the flexible 
boundary conditions involves signiflcantly less atoms than what is required by flxed bound- 
ary conditions methods. This highlights the applicability of flexible boundary conditions 
methods to modeling defect /interface interactions by DFT. 
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Appendix A: Evaluation of G ,j,(kx = 0) 
1. a = a' 

G luikx = 0) is obtained by taking the limit of Eqn. (IT2|) and Eqn. ( !T3|) as k^ — )■ 0: 

CL.(^- = 0)= / G--{Kz){e-''''^'''''^-^)dk. (Al) 



>(0) 
kf{0) q" 



+ / G''''ikJ)-^Mk,)dk, (A2) 

'fc,(0) «^z 



^2 

fe,(0) '^2 



-(/,(A;,) - l)dk, (A3) 



J. p^<^<^ / arctan(A;/(A;^)/A;a:) - arctan(A;j(A;^)/A;^) 7i_\ 

Note that since k = k^z, G {k) is evaluated along a constant /c-direction and therefore is a 

constant. The cut off function is 

f 1 < \kj < O.SA;^^ 

{ 12(1 - %\f - 16(1 - \k,\f 0.5A;^^^ < %\ < A;^^^ 

where fc™*" < Min(|A;j(0)| , A;/(0)) to insure that fc{kz) = at the Brillouin zone boundary. 
We isolate the k^ = point by dividing the integration path in Eqn. (lAip . Eqn. (lA2p and 
Eqn. flA3l) into three intervals 

[A;,(0), kf{0)] = MO), -e/2) U [-e/2, e/2] U (e/2, kj{0)] 

where e is sufficiently small. The first and third intervals do not contain the F-point and 
therefore their corresponding integrals are evaluated numerically without special treatments. 
To evaluate the integrals in Eqn. flAip and Eqn. flA2p over [—e/2, e/2], we use the small k^ 
leading order terms of G""{k) [14] and the exponential term 

e/2 

G"" {k,z){e-'^^^'''-''^'^ - l)dk, = 

-e/2~ 



L if + M^- M ^ ^-'^'^^ i"'^^'^ - ''^ - ''^^-^) ''^ 



G-%z,-zA-G'"^-^^^^]e 
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and 



.12 ^- 

-e/2 l^z 






e/2 ■ 



fcz 



G Jkl and G^^,{kz) are the elastic and discontinuity corrections and iG'^/kz appears only 



rrr 



in the case of a multiatom basis. G , and G""'^ are constants here[l4, 051. Also note that 

~ aa' ~aa i — -j 

fc{kz) = 1 over [—e/2, e/2]; hence the integral in Eqn. (1A3P equals zero over this interval. 

Taking e to be -^^^r — '— where A^div is the number of divisions in the discrete kz mesh 
we have 



^- ^ kfjO) - fc,(0) 

^al,a'l'y^) JVdiv 



1 1 



J2 [GZ'ikzz)e-''^^^'-^'^ 



G 



kl 



+ GL(0) 



+ G , 



aa' 



hiQ) kj{0)J- 



The first summation 



km - km 



N, 



div 



J2[G^c.a'{kzz)e~''^ 



i^l-^l') _ 



G 



kz^O 



i,2 



kf{0) km_ ^ (G^^^,(fc^5)(e-^^(^.-.) _!))+( G^"{kzz) - ^Ukz 



N. 



div 



kz¥=0 



kl 



G 



kl 



iUkz) - 1) 



is the numerical integration of all three integrals in Eqn. ( ]All) -( lA3l) over [A;j(0), — e/2) U 
(e/2, kf{Q)\. The last term G , ( pW — ^i^ ) is the evaluation of Eqn. flA4p . 



2. a^a' 



G ,{kx = 0) is obtained from the small kx expansion of Eqn. (Q and removing the k^ ^ 
term 



TT 



~aa'^ ' -T I ^ °^" 

t2 



Pll' 
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